Combine the salinity color plot with velocity vectors
%load_ext autoreload
%autoreload 2
import numpy as np
import pandas as pd
import xarray as xr
#
from schimpy import schism_mesh
# viz
import hvplot.pandas
import hvplot.xarray
import holoviews as hv
from holoviews import opts,dim
hv.extension('bokeh')
# widgets
import panel as pn
pn.extension()
# for advanced viz ops
import datashader
import holoviews.operation.datashader as hd
from holoviews.operation.datashader import datashade, rasterize
import warnings
warnings.filterwarnings('ignore')
varname = 'salinity'
dsv = xr.open_mfdataset('../tests/data/m1_hello_schism/outputs/out2d_*.nc', concat_dim='time', combine="nested",
data_vars='minimal', coords='minimal', compat='override')
dsv
<xarray.Dataset>
Dimensions: (time: 144, one: 1, nSCHISM_hgrid_node: 2639,
nSCHISM_hgrid_face: 4636,
nSCHISM_hgrid_edge: 7274,
nMaxSCHISM_hgrid_face_nodes: 4, two: 2)
Coordinates:
* time (time) datetime64[ns] 1999-12-31T16:20:00 ... 20...
SCHISM_hgrid_node_x (nSCHISM_hgrid_node) float64 dask.array<chunksize=(2639,), meta=np.ndarray>
SCHISM_hgrid_node_y (nSCHISM_hgrid_node) float64 dask.array<chunksize=(2639,), meta=np.ndarray>
SCHISM_hgrid_face_x (nSCHISM_hgrid_face) float64 dask.array<chunksize=(4636,), meta=np.ndarray>
SCHISM_hgrid_face_y (nSCHISM_hgrid_face) float64 dask.array<chunksize=(4636,), meta=np.ndarray>
SCHISM_hgrid_edge_x (nSCHISM_hgrid_edge) float64 dask.array<chunksize=(7274,), meta=np.ndarray>
SCHISM_hgrid_edge_y (nSCHISM_hgrid_edge) float64 dask.array<chunksize=(7274,), meta=np.ndarray>
Dimensions without coordinates: one, nSCHISM_hgrid_node, nSCHISM_hgrid_face,
nSCHISM_hgrid_edge,
nMaxSCHISM_hgrid_face_nodes, two
Data variables:
minimum_depth (one) float64 dask.array<chunksize=(1,), meta=np.ndarray>
SCHISM_hgrid (one) |S1 dask.array<chunksize=(1,), meta=np.ndarray>
depth (nSCHISM_hgrid_node) float32 dask.array<chunksize=(2639,), meta=np.ndarray>
bottom_index_node (nSCHISM_hgrid_node) int32 dask.array<chunksize=(2639,), meta=np.ndarray>
SCHISM_hgrid_face_nodes (nSCHISM_hgrid_face, nMaxSCHISM_hgrid_face_nodes) float64 dask.array<chunksize=(4636, 4), meta=np.ndarray>
SCHISM_hgrid_edge_nodes (nSCHISM_hgrid_edge, two) float64 dask.array<chunksize=(7274, 2), meta=np.ndarray>
dryFlagNode (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
elevation (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
depthAverageVelX (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
depthAverageVelY (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
dryFlagElement (time, nSCHISM_hgrid_face) float32 dask.array<chunksize=(72, 4636), meta=np.ndarray>
dryFlagSide (time, nSCHISM_hgrid_edge) float32 dask.array<chunksize=(72, 7274), meta=np.ndarray>ds = xr.open_mfdataset(f'../tests/data/m1_hello_schism/outputs/{varname}_*.nc', concat_dim='time', combine="nested",
data_vars='minimal', coords='minimal', compat='override')
ds
<xarray.Dataset>
Dimensions: (time: 144, nSCHISM_hgrid_node: 2639, nSCHISM_vgrid_layers: 11)
Coordinates:
* time (time) datetime64[ns] 1999-12-31T16:20:00 ... 2000-01-02T16:00:00
Dimensions without coordinates: nSCHISM_hgrid_node, nSCHISM_vgrid_layers
Data variables:
salinity (time, nSCHISM_hgrid_node, nSCHISM_vgrid_layers) float32 dask.array<chunksize=(72, 2639, 11), meta=np.ndarray>smesh = schism_mesh.read_mesh('../tests/data/m1_hello_schism/hgrid.gr3')
smesh.n_nodes(), smesh.n_edges(), smesh.n_elems()
(2639, 7274, 4636)
nodes = pd.DataFrame(smesh.nodes,columns=['x','y','z'])
nodes
nodes[varname] = ds[varname].values[0,:,0]
trimesh = hv.TriMesh((smesh.elems, hv.Points(nodes, vdims=varname)))
trimesh = trimesh.opts(
opts.TriMesh(cmap='rainbow4', cnorm='eq_hist', colorbar=True, node_alpha=0, edge_alpha=0, edge_color=varname, filled=True, height=400,
inspection_policy='edges', tools=['hover'], width=400))
trimesh.nodes.data
| x | y | index | z | salinity | |
|---|---|---|---|---|---|
| 0 | 56000.00 | -10400.00 | 0 | 0.500000 | 3.364312e-02 |
| 1 | 55831.58 | -10400.00 | 1 | 0.500000 | 1.095413e-02 |
| 2 | 56000.00 | -10350.00 | 2 | 0.944444 | 2.154961e-02 |
| 3 | 55663.16 | -10400.00 | 3 | 0.500000 | 2.462938e-08 |
| 4 | 55831.58 | -10348.07 | 4 | 1.023290 | 5.119039e-03 |
| ... | ... | ... | ... | ... | ... |
| 2634 | 55831.58 | 10348.07 | 2634 | 1.023290 | 5.119046e-03 |
| 2635 | 55663.16 | 10400.00 | 2635 | 0.500000 | 2.464225e-08 |
| 2636 | 56000.00 | 10350.00 | 2636 | 0.944444 | 2.154963e-02 |
| 2637 | 55831.58 | 10400.00 | 2637 | 0.500000 | 1.095410e-02 |
| 2638 | 56000.00 | 10400.00 | 2638 | 0.500000 | 3.364312e-02 |
2639 rows × 5 columns
vmag = np.sqrt(dsv.depthAverageVelX**2+dsv.depthAverageVelY**2)
vangle = np.arctan2(dsv.depthAverageVelY,dsv.depthAverageVelX)
vel = xr.Dataset({'mag':vmag,'angle': vangle})
vel## Calculate velocity vectors from depth averaged X and Y velocity vectors
<xarray.Dataset>
Dimensions: (time: 144, nSCHISM_hgrid_node: 2639)
Coordinates:
* time (time) datetime64[ns] 1999-12-31T16:20:00 ... 2000-0...
SCHISM_hgrid_node_x (nSCHISM_hgrid_node) float64 dask.array<chunksize=(2639,), meta=np.ndarray>
SCHISM_hgrid_node_y (nSCHISM_hgrid_node) float64 dask.array<chunksize=(2639,), meta=np.ndarray>
Dimensions without coordinates: nSCHISM_hgrid_node
Data variables:
mag (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
angle (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>def velocity_field(time, vector_size=1):
data = vel.isel(time=time)
vf = hv.VectorField((vel.coords['SCHISM_hgrid_node_x'],vel.coords['SCHISM_hgrid_node_y'],data.angle, data.mag))
vf.opts(opts.VectorField(pivot='tip', color=dim('Magnitude'), magnitude=dim('Magnitude').norm()*2000*vector_size, rescale_lengths=False))
return vf.opts(width=600)
velocity_field(70,1)
def salinity_mesh(time, level):
trimesh.nodes.data[varname] = ds[varname].isel(time=time, nSCHISM_vgrid_layers=level)
return rasterize(trimesh, precompute=False, dynamic=False,
aggregator=datashader.mean(varname)).opts(
colorbar=True, cmap='rainbow4', width=600, height=300, tools=['hover'])
salinity_mesh(0,0)
def update(time, level):
return salinity_mesh(time,level).opts(width=1000, height=500)*velocity_field(time, 1)
import panel as pn
pn.extension()
time_slider = pn.widgets.DiscreteSlider(name='Time',
options=dict(zip(ds.time.astype('str').values,
range(len(ds.time)))))
level_selector = pn.widgets.Select(name='nSCHISM_vgrid_layers',
options=dict(zip(ds.nSCHISM_vgrid_layers.astype('str').values,
range(len(ds.nSCHISM_vgrid_layers)))))
slider_row = pn.Row(time_slider, level_selector)
map_row = pn.Row(hv.DynamicMap(update,streams={'time':time_slider, 'level': level_selector}))
dash = pn.Column(slider_row, map_row)
dash.servable(f'HelloSCHISM: {varname}')#.show()